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METHOD OF FORMING A MODEL REPRESENTATIVE OF THE DISTRIBUTION 
OF A PHYSICAL QUANTITY IN AN UNDERGROUND ZONE, FREE OF THE 
EFFECT OF CORRELATED NOISES CONTAINED IN EXPLORATION DATA 

BACKGROUND OF THE INVENTION 

Field Of The Invention 

[0001] The present invention relates to a method of forming, from data obtained 
by exploration of a zone of a heterogeneous medium, a model representative of the 
distribution in the zone of a physical quantity (at least partly) free of the presence of 
correlated noises that may be contained in the data. 

[0002] The method applies, for example, to the quantification of the acoustic 
impedance in an underground zone. 


Description of the Prior Art 

[0003] The process of seeking a model that adjusts to experimental 

measurements has been developed in nearly all the fields of the sciences or 

i 

technology. Such an approach is known under various names: least-squares method 
for parameters estimation, for inverse problem solution. For a good presentation of 
this approach within the context of geosciences, one may for example refer to: 

Tarantola, A.: "Inverse Problem Theory: Method for Data Fitting and Model 
Parameter Estimation", Elsevier, Amsterdam, 1987. 

[0004] It can be noted that the term "least squares" refers to the square of the 
norm in the data space for quantifying the difference between the response of a 
model (which is the image of the model by a previously selected modelling operator) 
and the data using, a cost function that has to be minimized to solve the problem. 
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Using the square of the norm to define the cost function is just a practical 
convenience, and it is not fundamentally essential. Besides, many authors use, for 
various reasons, another definition of the cost function but this definition remains 
based on the use of the norm, or of a semi-norm, in the data space. Finally, 
considerable latitude exists for selecting the norm (or the semi-norm) in the data 
space (the use of the Euclidean norm is not required). In the case of noise- 
containing data, the solution can substantially depend on the choice made at this 
stage. More developments concerning this problem, are referred to in the following 
publications: 

Tarantola, A.: "Inverse Problem Theory: Method for Data Fitting and Model 
Parameter Estimation", Elsevier, Amsterdam, 1987; Renard and Lailly, 2001; 
Scales and Gersztenkorn, 1988 ; Al-Chalabu, 1992. 

[0005] The measuring results often contain errors. Modelling noises add further to 
these measuring errors when the experimentation is compared with modelling 
results: modellings are never perfect and therefore always correspond to a simplified 
view of reality. The noise will therefore be described hereafter as the superposition 
of: 

non correlated components (white noise for example), 

correlated components, which means that the existence of a noise on a 
measuring sample translates into the existence of a noise of equal nature on certain 
neighbouring measuring points ; modelling noises typically belong to this category. 

[0006] When the data contain correlated noises, the quality of the model 
estimated by solving the inverse problem can be seriously affected thereby. As 
already mentioned, no modelling operator is perfect. It is therefore the work of the 
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whole community of the people involved in the identification of parameters describing 
a model which is hindered by the existence of correlated noises. Among these 
people, seismic exploration practitioners are among the most concerned ones: in 
fact, their data have a poor or even very bad signal-to-noise ratio. This is the reason 
why correlated noise filtering techniques are an important part of seismic data 
processing softwares. The most conventional techniques use a transform (Fourier 
transform for example) where signal and noise are located in different areas of 
space, thus allowing separation of the signal from the noise. A general presentation 
of the conventional methods for noise elimination on seismic data is discussed in the 
book by Yilmaz: 

Yilmaz, O. 1987: "Seismic Data Processing", Investigation in Geophysics 
No.2, Society of Exploration Geophysicists, Tulsa, 1987. 

[0007] However, these techniques are not perfect: they presuppose that a 
transformation allowing complete separation of the signal and of the noise has been 
found. In this context, correlated noises are particularly bothersome because they 
can be difficult to separate from the signal (which is also correlated) and can have 
high amplitudes. It is therefore often difficult to manage the following compromise: 
either the signal is preserved, but a large noise residue remains, or the noise is 
eliminated, but then the signal is distorted. These filtering techniques can be 
implemented before solving the inverse problem: they constitute then a 
preprocessing applied to the data. The quality of the inverse problem solution then 
widely depends on the ability of the filters to eliminate the noise without distorting the 
signal. 

[0008] This approach is introduced by Nemeth et al. in the following publication: 
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Nemeth T. et al. (1999), "Least-Square Migration of Incomplete Reflection 
Data", Geophysics, 64, 208-221 which constitutes an important advance for the 
inversion of data disturbed by high-amplitude correlated noises. The authors propose 
eliminating the noise by solving an inverse problem: after defining the space of the 
correlated noises as the image space of a vector space B (space of the noise- 
generating functions) by a linear operator T and the signal as the image space of a 
vector space M (models space) by a modelling operator F (assumed to be linear), 
they seek the signal in F(M) and the noise in T(B) whose sum is as close as possible 
(in the sense of the Euclidean norm or, on a continuous basis, of the norm L 2 ) to the 
measured datum. This technique constitutes an important advance insofar as it 
allows elimination (or at least very substantial reduction) of high-amplitude correlated 
noises (surface waves for example) that are difficult to separate from the signal by 
means of conventional filtering techniques. However, according to the authors, an 
increase by an order of magnitude of the computing time required for solution of the 
conventional inverse problem, that is without seeking the correlated component of 
the noise, is the price to pay for these performances. Furthermore, the result 
obtained with the method is extremely sensitive to any inaccuracy introduced upon 
definition of operator T: this is the ineluctable compensation for the high aptitude of 
the method to discriminate signal and noise. 


SUMMARY OF THE INVENTION 

[0009] The method according to the invention allows to estimate, from data 
obtained by exploration of an underground zone, a model representative of the 
distribution, in the zone, of at least one physical quantity, this model being 


4 


' 'SUBSTITUTE SPECIFICATION 612.42948X00 

substantially free of the presence of correlated noises that may be contained in the 
data. The method essentially comprises the following stages: 

a) acquisition of measurements giving information about certain physical 
characteristics of the zone by following a predetermined experimental protocol, 

b) specification of a modelling operator which associates, with a model of 
each physical quantity, synthetic data that constitute the response of the model, the 
measurements and the synthetic data belonging to a data space, 

c) for each correlated noise referenced by a subscript j ranging from 1 to 
J, selection of a modelling operator which associates a correlated noise with a noise- 
generating function belonging to a predetermined space of noise-generating 
functions, 

d) specification of a norm or of a semi-norm in the data space, 

e) specification of a semi-norm in the space of the noise-generating 
functions for which each noise modelling operator establishes substantially an 
isometric relation between the space of the noise-generating functions space and the 
data space, 

f) definition of a cost function quantifying the difference between the 
measurements on the one hand and the superposition of the model response and of 
the correlated noises associated with the noise-generating functions on the other 
hand, and 

g) adjustment of the model and of the noise-generating functions by 
minimizing the cost function, by means of an algorithmic method taking advantage of 
the isometry properties of the noise modelling operators. 


5 


* Substitute specification 612.42948x00 

[0010] According to a first implementation mode, the distribution as a function of 
the depth of the acoustic impedance in the medium is sought, the correlated noises 
affecting the data are tube waves identified each by parameters characterizing their 
propagation, the measured data are VSP type data obtained by means of pickups 
suited to detect the displacement of particles in the medium in response to a 
localized seismic excitation, the location of the pickups, the recording time and the 
time sampling points being defined, and the modelling operator selected associates 
the synthetic data with an acoustic impedance distribution as a function of the 
evaluated depth in traveltime and with the vertical stress measured as a function of 
time at the depth of the first pickup. 

[0011] The cost function selected to quantify the difference between the 
measurements on the one hand and the superposition of the model response and of 
the correlated noises associated with the noise-generating functions on the other 
hand is for example the square of the semi-norm of this difference in the data space. 

[0012] According to an embodiment, adjustment of the model and of the noise- 
generating functions is for example obtained by means of a relaxation method for 
eliminating the unknowns corresponding to each correlated noise generating 
function, this relaxation method being implemented within the iterations of a quasi- 
Newtonian algorithm for calculation of the model. 

[0013] According to an embodiment, numerical calculation of the image of a 
model by this operator is carried out for example by numerical solution of the 1D 
waves equation for the model considered, by selecting values taken by the 
displacement of the particles at the locations of pickups and at the previously defined 
time sampling points, and by applying an operator likely to compensate for the 
spherical divergence and attenuation effects. 
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[0014] According to an embodiment, the numerical noise modelling operator is a 
finite-difference centered numerical scheme for discretizing the noise transport 
equation, and the noise-generating function involved as the initial condition along the 
edge of the zone associates with each correlated noise a space of noise-generating 
functions comprising support time functions in a given time interval. 

[0015] Various examples of norms or semi-norms selected for the data space and 

the space of the noise-generating functions are given in the detailed description 

1 

hereunder. 

[0016] According to a second implementation mode, the distribution of 
disturbances, in relation to a previously selected reference medium, of the 
impedance and of the velocity in the zone of the medium is sought, the correlated 
noises affecting the data are due to multiple reflections whose kinematics and 
amplitude variations with the offset have been previously estimated, the measured 
data are picked up by seismic surface pickups, the location of the pickups, the 
seismic excitation mode, the recording time and the time sampling points being 
defined, and the modelling operator being defined via a linearization of the waves 
equation around the reference medium. 

[0017] The cost function selected to quantify the difference between the 
measurements on the one hand and the superposition of the model response and of 
the correlated noises associated with the noise-generating functions on the other 
hand is for example the square of the semi-norm of this difference in the data space. 

[0018] According to an embodiment, adjustment of the model and of the noise- 
generating functions is obtained by means of a block relaxation method for 
eliminating the unknowns corresponding to each correlated noise generating 
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function, this relaxation method being implemented within the iterations of 
conjugate gradient algorithm for calculation of the model. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

[0019] Other features and advantages of the method according to the invention 
will be clear from reading the description hereafter of an embodiment given by way 
of non limitative example, with reference to the accompanying drawings wherein: 

Fig. 1 shows an example of data obtained by means of a vertical seismic 
prospecting method (VSP), 

Fig. 2 shows a model of distribution of the acoustic impedance with depth, 

Fig. 3 shows synthetic VSP data obtained on the basis of the impedance 
model of Fig.2, 

Fig. 4 shows an example of VSP data contaminated by a single correlated 

noise, 

Fig. 5 shows an example of VSP data contaminated by two correlated noises, 

Fig. 6 shows the impedance distribution as a function of depth, obtained by 
inversion of the noise-containing data of Fig.4, 

Fig. 7 shows the inversion residues (differences between the data of Fig.4 and 
the seismic response of the model of Fig.6), 

Fig. 8 shows the impedance distribution as a function of depth, obtained by 
inversion of the noise-containing data of Fig.5, 

Fig. 9 shows the corresponding inversion residues (differences between the 
data of Fig.5 and the seismic response of the model of Fig. 8), 

Fig. 10 shows the impedance distribution as a function of depth, obtained by 
inversion of the noise-containing data of Fig.4 by seeking the correlated noise in a 
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form of the superposition of two correlated noises having inaccurate propagation 
properties, that is different from those of the noise appearing in the data of Fig.4, 

Fig. 11 shows the corresponding inversion residues (differences between the 
data of Fig.4 and the seismic response of the model of Fig. 10), 

Figs. 12A and 12B respectively show seismic data with multiple reflections 
and the seismic response of the model obtained after conventional linearized 
inversion, 

Figs. 13A and 13B respectively show an example of impedance distribution 
and of propagation velocity for which the seismic data of Fig.12A constitute the 
seismic response, and 

Figs. 14A and 14B respectively show the comparison of the seismic 
responses of the models obtained by conventional inversion (Fig. 14A) and by 
inversion according to the proposed method involving picking of the moveout of the 
multiple and an estimation of the amplitude variations with the offset (Fig. 14B). 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
[0020] Formulation of the problem 

[0021] Data is available resulting from a sampled measurement of a function. 
This function depends on several variables (space or space-time variables for 
example). 

[0022] These data, called d, correspond to measurements carried out to obtain 
information about a model m. They contain various noises: 

correlated noises (or a superposition of correlated noises), 
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non correlated noises. 
[0023] The problem is to determine quantitatively model m (or functions of this 
model) from data d. 

[0024] A modelling operator F is selected (linear or not) which associates with 
model m the response of the model This operator actually models the real physical 
phenomenon only imperfectly. This is essentially (but not exclusively) the reason 
why correlated noises appear in the data: these noises correspond to a signal 
related to the model but the relation between them appears to be too complex to be 
included in modelling operator F, which should keep relatively simple for various 
reasons. 

[0025] Under such conditions, the noise is modeled by introducing one or more 
noise-generating spaces Bj (j ranging from 1 to J) and the associated noise 
modelling operators called Tj. In order to find model m from data d, seeking this 
model is proposed as the solution to the optimization problem. 

min j C(m, p l , p 2 y ...fij ) = If(iw) + YJfi} - J (1) 

where || jj^ is a norm (or possibly a semi-norm) in the data space. 

[0026] The above relationship differs from the approach proposed by Nemeth et 
al. (1999) in that: 

modelling operator F is not necessarily linear, 

considering a superposition of correlated noises of different types is justified 
(in the sense that they are modelled by different operators referenced by subscripts 
j), at the price of a marked increase in the number of calculations concerning the 
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number of unknowns (which particularly complicates the solution of the optimization 
problem) as well as the number of noise modellings to be carried out, 

the possibility of selecting norm j| \\ D or possibly the semi-norm to suit our 
convenience in the data space is reserved. 

[0027] The choice that is made for |j \\ D will be imposed by considerations linked 

with the efficiency of the solution method (described below), this efficiency allowing 
dealing with the case of a correlated noise superposition. This is an important 
opening: besides the possibility of processing data containing correlated noises with 
complex characteristics, in this way the difficulties are overcome which are 
encountered in Nemeth et al.'s approach, linked with the high sensitivity of the result 
to an inaccuracy introduced upon definition of the noise modelling operator. To 
overcome this difficulty, the noise is sought by superposition of noises associated 
with approximate modelling operators Tj. 

[0028] It can be noted that selection of the square of norm |j || D involved in the 

optimization problem is not essential: the solution is not changed by selecting 
instead a function which, like the square function, is an increasing function on 9?* (or 
a decreasing function provided that the min is changed into the max). 

[0029] The solution method 

[0030] The method suitably selects: 

norm (or semi-norm) || in the data space, 

norms || \\ D in each space Bj 
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so as to obtain a high-performance algorithmic solution to optimization problem (1). 
The expression "suitably select" means see that each operator Tj constitutes an 
isometry (or an isometry approximation) for norms | | and || fl^ respectively in the 

initial space and in the end space. The main idea in the determination of norms, 
allows making each operator Tj to be isometric, is based on the existence of energy 
conservation type equalities verified by the solutions to the partial differential 
equations (or of discrete energy equalities for the solutions to certain numerical 
schemes used for discretization of these equations). 

[0031] If such choices are made, it is possible to make the cost function 
minimization very simple by eliminating the unknowns associated with the correlated 
noise generating functions (determination of Schur's complement). This elimination 
is peformed by means of a suitable algorithm such as the block relaxation method, a 
block being associated with a noise-generating function, or the conjugate gradient 
method with symmetrical-block Gauss-Seidel preconditioning, etc., all of them 
algorithmic implementation methods well-known in the art, presented for example 
by: 

Golub, G.H. et al. (1983) "Resolution Num6rique des Grands Systems 
Lin6aires" (Eyrolles). 

[0032] Simplified minimization of the cost function allows considerable reduction 
of the time required for numerical solution. 

[0033] First implementation example applied to 1D inversion of VSP data 
contaminated bv a tube wave 
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[0034] VSP (Vertical Seismic Profile) data belong to conventional well survey 
acquisitions. VSP data have multiple applications, one of the main ones being their 
use for establishing a connection between seismic surface data and logging 
measurements. The inversion of VSP data as described for example by Mac6, D. 
and Lailly, P. (1984) A Solution of an Inverse Problem with 1D Wave Equation 
Applied to the Inversion of Vertical Seismic Profile ; Proc. of the 6 th Intern. Conf. In 
"Analysis and Optimisation of Systems", Nice, 2, 309-323 has been developed to 
that end. In the'1D inverse problem, the subsurface model is assumed to be laterally 
invariant and that a plane-wave excitation propagates vertically. The unknowns of 
the problem are the acoustic impedance distribution as a function of depth 
(measured in vertical traveltime and over an interval starting at the depth of the first 
receiver) and the seismic excitation mode (precisely the time function characterizing 
the boundary condition at the depth of the first pickup). This is a non linear inverse 
problem, the VSP data depending non linearly on the impedance distribution. 

[0035] VSP data are often contaminated by a correlated noise of very high 
amplitudefrom a fluid wave. This wave is guided in the mud that invades the well. Its 
propagation velocity is low in relation to the propagation velocity in the rocks 
surrounding the well. Furthermore, the propagation velocity of the fluid wave 
depends on the frequency (dispersive propagation). Finally, the fluid wave reflects 
notably at the well bottom and furthermore gives rise to another propagation mode. 
Typical VSP data are given in the figure below: the main characteristics of the tube 
wave can be seen. Its amplitude and the extension of the contaminated zone make it 
difficult to use the signal contained in the many seismic samples contaminated by 
the fluid wave. 
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[0036] Hereafter the implementation of the method for 1 D inversion of VSP data 
is described. 

[0037] Experimental protocol 
[0038] It specifies: 

the various depths at which the receivers are positioned in the well (assumed 
to be vertical): in experiments, the pickups cover the range of depths measured in 
one-way time [x m i n = 0.05s ; Xmax = 0.2s], with a receiver every Ax=1 ms (in one-way 
time) ; thus samples numbered from 0 to 1 are obtained, 

the physical nature of the measurement performed by the pickups: in 
experiments, the pickups measure the vertical displacement resulting from the wave 
propagation, 

the recording time: in experiments, the pickups measure the vibrational state 
over the time interval [0, T=1.5s], these data being sampled every At=1 ms. Thus 
samples numbered from 0 to N are obtained. 

[0039] The data di n is the data sample recorded at the depth Xmin + i Ax and at the 
time n At. Of course: T = N At and x^ax = Xm in + I Ax. 

[0040] Acquisition of data in accordance with the experimental protocol 

[0041] Experiments were carried out from synthetic data, that is obtained by 
numerical modelling, itself carried out in two stages: 

[0042] Creation of synthetic data containing no noise: this creation was performed 
by numerical solution of the 1D acoustic wave equation, the model being the 
acoustic impedance distribution shown in Fig,2 (the depth is measured in one-way 
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time) and the excitation mode (Neumann's boundary condition at the depth 0) being 
a conventional Ricker wavelet. 

[0043] The synthetic VSP data resulting from this modelling are shown in Fig. 3: 
the vertical axis representing the observation time is graduated from 0 to 1 .5 s and 
the horizontal axis representing the depth (in one-way traveltime) of the various 
pickups is graduated from 0.05 to 0.2 s. 

[0044] Addition of noise to the data 

[0045] The VSP data are disturbed thus calculated by adding thereto, on the one 
hand, a random noise (but filtered in the seismic frequency band) of rather high 
amplitude and, on the other hand, one or more correlated noises of very high 
amplitude. Each noise will propagate downwards with a low velocity which depends 
on the frequency (dispersive propagation): these correlated noises are supposed to 
represent fluid waves. They have been modelled by numerical solution (use of the 
finite-difference centered scheme) of a transport equation with a constant 
propagation velocity. Dispersion of the waves has been obtained using large 
discretization intervals in the finite-difference scheme. In the case of the 
superposition of several correlated noises, the propagation velocities associated with 
each noise are different. Figures 4and 5 show the VSP data contaminated by a 
single correlated noise (Fig.4) and by the superposition of two correlated noises 
(Fig.5). 

[0046] Specification of the modelling operator 

[0047] Modelling is carried out by solving the wave equation: 

<r(x) £f - ± (cr(x) = 0 in ^ ,+co[ x [0, T] 
ot ox ox 
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with Neumann's boundary condition: 

ax 

and the initial conditions: 

[0048] Definition of the modelling operator first requires definition of the 
observation operator which formalizes the experimental protocol described in §1. 
The observation operator thus is the operator which associates with function y(x,t), 
solution to the above system, samples y(X|,t n ) for Xj=Xmin + i Ax and t n =n At, i=0,...,l ; 
n=0,...,N. 

[0049] The modelling operator identified as F(m) thus is the (non linear) operator 
which associates with function pair m=(a(x),h(t)) samples y(x,,t n ) for X|=Xmin + i Ax and 
t n =nAt, i=0 l;n=0,...,N. 

[0050] Specification of the modelling operator associated with each correlated 
noise 

[0051] The procedure is specified herein that can be used for modelling each 
correlated noise. Each correlated noise will be characterized by correlation directions 
specific thereto, which are assumed to be known, at least approximately: these 
correlation directions are specified by means of a field of correlation vectors dj(x 9 i) 

of components (q x (x,t), q'(x,t)) which has to be defined at any point of the domain 
[Xmin, Xmax]X[0,T]. In an equivalent way, the correlation directions can be specified by 
means of correlation lines that will represent the field lines associated with Cj(x 9 t). 
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Such correlation lines can be obtained according to the technique described in US 
Patent 4,972,383 filed by the applicant, from a pick of several phases characteristic 
of the noise, the information being extended to the whole domain [Xmm, Xmax]X[0,TJ by 
means of conventional interpolation and/or extrapolation procedures. Specification of 
the correlation directions amounts to specifying the propagation velocity distribution 
of the noise q(x,t) that is related to the correlation vector by the relation 


variations during propagation of the wave (a different situation is presented for the 
2 nd application example). Modelling of a correlated noise as specified above is based 
on the following observation: a wave propagating along correlation lines without 
amplitude change is solution to the transport equation: 


[0052] In this context, introduced are: 

a space Bj of noise-generating functions which will be the space of support 
functions ft(t) in [tj mln ? tj max ] c [0,T] and that is extended by 0 in the whole interval 

[o/n. 

noise modelling operator Ty. 

Tj.fijiOeBj^bjix^eD 

solution to the transport equation 
Vbj(x 9 02j{x 9 t) - 0 in L,^]x[0,r] 

and meeting the initial conditions: 



assumed that there are no wave amplitude 
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'*>.(*,/ = 0) = 0 

[0053] It can be noted that the geometry of the correlation lines (or, which 
amounts to the same thing, the propagation velocity distribution) cannot be any 
geometry: the correlation lines cannot intersect, or the transport equation solution 
problem would not be correctly set. For the same reasons, the correlation lines 
cannot intercept twice the edges on which the initial conditions are specified. The 
latter consideration can lead (in order to model a correlated noise corresponding to 
an upgoing wave) to specify the initial condition no longer at x=Xniin but at x=Xma X or 
even at an intermediate value, even if it means decomposing Cauchys problem into 
two problems associated with the domains located on either side of this intermediate 
edge. 

[0054] It is furthermore assumed that the geometry of the correlation lines is such 
that (more general situations can be considered, but the following results then take a 
different form): 

oj(x y t) = £j(x) x r(t) 
ou encore : <fj{x,t) = £ s (x) et c^(x f t) = ^ 

[0055] Finally, it is assumed that the pencil of correlation lines that intercept at 
x=x min interval [tj mln ,tj max ] do not intercept at t=T interval [x^in^aj. 

[0056] In fact, the transport equation is solved by means of a numerical scheme. 
It is for example possible to use the conventional finite-difference centered scheme. 
Therefore discretization intervals AXj' and Atj are selected (which are not necessarily 
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equal to Ax and At, but which are assumed to be submultiples of these quantities, 
even if more general situations can be considered) and a grid is introduced whose 
nodes are the points of coordinates (x^in + i'Axj. n'At/) with i' and n' e N. The 
conventional finite-difference centered scheme explained hereunder allows, from the 
initial conditions, to calculate step by step the different values taken by function 
bj(x,t) at the various nodes of the grid (to simplify the notations, subscript j 
associated with the correlated noise considered is omitted). 


[0057] In the above formula, quantity txr n represents the evaluation of quantity a 
(a is a generic term) at the point of coordinates (Xm in + Tax/, n'At/). 

[0058] Of course, the use of discretization intervals Ax/ and Atj sufficiently small in 
relation to the wavelength is required if the numerical scheme is to give a precise 
approximation of the function solution to the transport equation. 

[0059] However, the use of larger discretization intervals (but still smaller than the 
wavelength) allows modelling more complex propagation phenomena by making the 
wave propagation velocity dependent on the frequency: dispersive propagations 
such as the propagation of tube waves can thus be modelled. Selection of suitable 
discretization intervals for modelling a particular mode of the tube wave can be 
made: 



1 



.n' + L 


2Az' r 
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[0060] By determining, after examining the data, the propagation velocity of this 
mode (which depends here not only on space and time, but also on the frequency) ; 
tomography techniques can be used to that end, such as those presented in the 
following publication: 

Ernst, F. et al. (2000) ; Tomography of Dispersive Media ; J. Acoustic Soc. 
Am., 108, 105-116. 

[0061] By knowing the dispersive propagation properties associated with the 
numerical scheme, properties that follow from an analysis of the numerical 
dispersion relation, as described for example in the following publication: 

Alford, R.M. et al. (1974), Accuracy of Finite Difference Modeling of the 
Acoustic Wave Equation ; Geophysics, 6, 834-842. 

[0062] Specification of a norm or of a semi-norm in the data space 

[0063] A first choice takes as the norm in the data space any discrete norm which 
is meant to be an approximation (via a quadrature formula) to the norm defined on a 
continuous basis by: 

i 

[0064] The semi-norm, which is a better choice as discussed in the next 
paragraph may also be used: 
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where u is any vector of the data space (i and n are the indices representing the 
sample numbers respectively in space and in time). 

[0065] Besides, other choices are also possible. 

[0066] Specification, in each space Bj of noise-generating functions, of a norm for 
which each noise modelling operator 7} establishes an isometric relation (or the 
approximation of an isometric relation) between the noise-generating functions 
space and the data space (provided with the semi-norm as defined in §5) 

[0067] The procedure herein specifies defining the norm in the noise-generating 
functions space associated with each correlated noise. This procedure being the 
same for each noise, the description is generic and leaves out, in the formulas, 
subscript j associated with the noise considered. 

[0068] A first choice provides space of each noise-generating functions B with 
any discrete norm meant to be an approximation (via a quadrature formula) to the 
norm defined on a continuous basis by: 



in which case, using the fact that, by means of the hypotheses set out in §4, 
bj(x,T)=0 for x e [Xmin.XmaxL linear operator Tj performs an isometry between Bj and 
D. However, since calculation of Tjft is carried out numerically, the isometric 
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character of operator T will not be perfect but only approximate, the approximation 
being all the better as discretization intervals Ax, At, Ax/, At)' are small. The fact that 
only approximately an isometry is available leads to lower performances concerning 
the algorithmic implementation presented hereafter. 

[0069] A better choice provides each noise-generating functional space Bj with 
the semi-norm: 


[0070] In this case, operator Tj rigorously performs an isometry between Bj and D, 

at least if u f N = OVi = 0,1-1 . This results from the discrete energy equality: 

[0071] an equality valid in the case (simplified for presentation's sake) where Ax' 
= Ax and At 1 = At and if £. 1/2 = £1/2 has been defined. 

[0072] Definition of the cost function 

[0073] The cost function is defined by the formula as follows: 
C(«,>? 1 ,/? 2i ...A) = |F(w) + 2Vy- rf 
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[0074] Seeking the model and the noise-generating functions minimizing the cost 
function, this search being carried out by an algorithmic method taking advantage, 
explicitly or implicitly, of the isometry properties of the noise modelling operators (see 
6) 

[0075] The algorithmic method uses the specificity of Schur's complement, a 
specificity associated with the isometric character of operators Tj. This can be done 
for example by realizing that minimizing C(m,Pi, p 2 ,..., pj) amounts to minimizing 
C\m)^C(m 9 fi x (mlfi 2 (m) y ^fij(m)) where (fi l (m) 9 fi 2 (m) t ...fi J (m)) minimizes C(m,p 1f 

p 2 Pj) for given m. It can be noted that the Hessian associated with this quadratic 

form is Schur*s complement. C'(m) can be minimized by means of any optimization 
method such as the BFGS version of a quasi-Newtonian method. 

[0076] Determination of (A(m),^ 2 (/w),...^ y (m) is very easy as a result of the 
isometric character of operators T j( notably if choices (2) and (3) have been 
respectively selected to define \\\\ D and \\\\ g . If there is only one type of noise (J=1) 

and if the isometry is perfect, determination of p x {m) simply requires evaluation of 
the gradient of the quadratic form. If there is a superposition of several noises (J>1 ), 
various algorithms can be considered to take advantage of the isometric character of 
operators Tj. (j0,(m),jS 2 (m),...£ y (ro) can for example be determined by means of a 
block relaxation method (Gauss-Seidel alias), a block being associated with a group 
of unknowns /? y (m), a relaxation iteration simply requiring (still in the case of a 

perfect isometry) a single evaluation of the gradient of the quadratic form associated 
with the block considered. For problems made difficult by the proximity of the 
correlation directions of the various noise types, a preconditioned conjugate gradient 
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method can be considered, the preconditioning matrix being of symmetrical block 
(relaxation alias) Gauss-Seidel type. Then again, preconditioning is very easily 
implemented and solution of the problem associated with a block only requiring 
evaluation of the gradient of the quadratic form associated with the block considered. 

[0077] Figs. 6 to 11 show the results obtained with three types of experiments 
using the method. 

[0078] In the case of an inversion of noise-containing data (Fig.4) both by 
superposition of a correlated noise and of a random noise, Fig.6 gives the 
impedance distribution obtained, and Fig.7 shows the inversion residues. 

[0079] In the case of an inversion of noise-containing data (Fig.5) both by 
superposition of two correlated noises and of a random noise, Fig.8 gives the 
impedance distribution obtained, and Fig.9 shows the inversion residues. 

[0080] In the case of an inversion of noise-containing data (Fig.4) both by 
superposition of a correlated noise and of a random noise where the propagation 
properties of the fluid wave are not precisely known, it is natural to seek a correlated 
noise model which superposes two correlated noises having inaccurate propagation 
velocities but by which the true propagation velocity of the fluid wave is framed. 
Fig. 10 gives the impedance distribution obtained and Fig. 11 shows the inversion 
residues. 

[0081] Second example of application of the method to determination of a model 
by linearized inversion of multi-offset data contaminated by multiple reflections 

[0082] This application differs from the first one mainly in that it illustrates the 
possibilities afforded by the method for taking into account noise amplitude variations 
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along the correlation directions. Another difference lies in the propagation properties 
of the correlated noises, propagation of the multiple reflections being not dispersive. 

[0083] Linearized inversion is one of the sophisticated methods used by 
geophysicists to obtain a quantitative model of subsurface heterogeneities, this type 
of information being important notably for reservoir characterization. The data given 
here are surface seismic data. This method is described for example in the following 
publication: 

Bourgeois, A. et al. (1989) The Linearized Seismic Inverse Problem: An 
Attractive Method for a Sharp Description of Strategic Traps": 59 th Ann. Intern. Mtg. 
Soc. Expl. Geoph. Expanded Abstract, pp.973-976. 

[0084] Furthermore, it is necessary to have a reference model, consisting (in the 
context of an acoustic model) of a velocity distribution and of an acoustic impedance 
distribution: it is the model around which a linearization will be carried out in the 
waves equation (Born's approximation). In order to simplify the presentation, the 1D 
case is taken but this hypothesis is all the less essential as the practical applications 
mainly concern 3D seismic methods. In the 1 D context, all of the seismic information 
is contained in a single shotpoint. Besides, the reference model is also one- 
dimensional: it only depends on the depth. The unknowns of the problem are the 
acoustic impedance disturbance and velocity distribution as a function of depth. On 
the other hand, unlike the first application example, the seismic excitation mode (and 
notably the seismic wavelet) are assumed to be known. It is an inverse problem 
made linear by means of linearization. 

[0085] Surface seismic data are often contaminated by multiple reflections which 
may have a rather high amplitude. Their kinematics, which is characterized by 
moveout in the 1D context, is generally quite different from that of primary reflexions: 

26 


' Substitute specification 61 2.42948x00 

these two types of waves have generally travelled through geologic layers having 
different propagation velocities. This is the reason why linearized inversion is in itself 
an efficient technique for attenuating multiple reflections: the kinematics of primary 
reflections (the events modelled by Born s approximation) is entirely defined by the 
velocity distribution in the reference model: the modelled events therefore cannot 
adjust to multiple reflections if the moveout thereof differs from that of the primary 
reflections. However, because of the stationarities at the zero offset, there is in 
practice a partial adjustment and the linearized inversion results appear to be 
contaminated by the presence of the correlated noise consisting of the multiple 
reflections. The seismic data of Fig.12A have been modelled by taking account of 
the multiple reflections (there are only 3 primary reflections which reach the zero 
offset successively at the times: 500; 1000; 1650 ms, the 2 nd arrival having a 
particularly high amplitude). 

[0086] In the seismic response of the model obtained after conventional 
linearized inversion (see Fig.12B), the high-amplitude multiple arriving at 
approximately 1500 ms at the zero offset has been only weakly attenuated and the 
model is contaminated by the multiple. It is then natural to try to use the possibilities 
afforded by Nemeth's method in order to better discriminate multiples and primaries, 
notably with the improvements provided by our method. 

[0087] Hereafter the implementation of the method for 1 D linearized inversion of 
surface seismic data is described. 

[0088] Experimental protocol 
[0089] It specifies: 
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the seismic excitation mode: in experiments, the seismic source is a source 
point time-modulated by a function (seismic wavelet) denoted by w(t), 

the various locations where the receivers are positioned: in experiments, the 
pickups are arranged at the depth 10m and cover the offset interval [Xmi n =0, 
x ma x=1500ml, with a receiver every Ax=100m ; thus obtained are seismic traces 
numbered from 0 to I, 

the pickups measure the pressure field resulting from the wave propagation, 

the recording time:: in experiments, the pickups measure the vibrational state 
over the time interval [0, T=1 800ms], these data being sampled every At=1 ms. Thus 
samples are obtained numbered from 0 to N. 

[0090] The quantity d* n us referred to as the data sample recorded for the offset 
x min + i Ax and at the time n At. The relationship: T = N At and x max = Xmin + I Ax 
exists. 

[0091] Acquisition of data in accordance with this experimental protocol 

[0092] Experiments were carried out from synthetic data: these data were 

JL^ _ v - VP - 8(x, z - z s )w(t) in 91 x 91 ** x[0,j] 
oc dt c 

obtained by numerical resolution (finite-difference method) of the 2D wave equation: 

with the boundary condition: 
P(z=0)=0 

and the initial conditions: 
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/>(< = 0) = ^(/ = 0) = 0 
at 

an equation wherein: 

x, z, t respectively designate the lateral coordinate, the depth and the time, 

a and c respectively represent the acoustic impedance and propagation 
velocity distributions (functions of the depth). 

[0093] The wavelet selected is a conventional Ricker wavelet of central frequency 
30 Hz and the impedance a and velocity c distributions are given in Figs.13A, 13B. 
The seismic response of this model therefore is the data represented in Fig.12A. The 
quantity F N l is referred to as the operator which associates the seismic data with 
function pair (a(z).c(z)). 

[0094] Specification of the modelling operator 

[0095] A reference model is selected which is described by function pair 
(a 0 (z),c 0 (z)). In the experiments presented hereafter, co(z)=c(z) and a 0 (z)=cte (=1) 
have been selected. The Jacobian operator of operator Fnl is chosen as the 
modelling operator at point (a 0 (z),Co(z)). In fact, since Cq(z)=c(z) has been selected, 
the only unknown is 8a(z)=a(£)-a 0 (z) and only the corresponding component of the 
Jacobian is important. 

[0096] The modelling operator, which is referred to as F(8m), thus is the (linear) 
operator which associates with function pair 5m=(8a(z),8c(z)) the samples 8y(x i ,t n ) 
which are the components of vector F (8m) for x i =x m i n + i Ax and t n =n At, i=0,...,l ; 
n=0,...,N. 
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[0097] Specification of the modelling operator associated with each correlated 
noise 

[0098] The procedure is specified herein that can be used for modelling each 
correlated noise. Each correlated noise is characterized by correlation directions 
specific thereto, which are assumed to be known, at least approximately: these 
correlation directions are specified by means of a field of correlation vectors Cj(x 9 t) 

of components (Cj x (x,t), c/(x,t)) which has to be defined at any point of the domain 
[Xmin. XmaJXfO.T]. In an equivalent way, the correlation directions can be specified by 
means of correlation lines that represent the field lines associated with c y Such 

correlation lines can be obtained according to the technique described in US Patent 
4,972,383 filed by the assignee, from a pick of several phases characteristic of the 
noise, the information being extended to the whole domain [Xmj n , x max ]X[0,T] by 
means of conventional interpolation and/or extrapolation procedures, as described in 
the aforementioned patent. For this second application example wherein elimination 
of the aforementioned multiple is desired, this multiple is picked (the amplitude peak 
for example), which enables defining the arrival time variations as a function of the 
offset and defines a correlation line continuum by simple vertical translation of the 
chosen line (under these conditions, vector c(x,t) does not depend on t). 
Specification of the correlation directions amounts to specifying the propagation 
velocity distribution of the noise Cj(x,t) that is related to the correlation vector by the 


[0099] Unlike the 1 st application example described above, here noise amplitude 
variations are along the correlation line. The assumption is that these amplitude 
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variations can be estimated from a measurement on the data or by means of 
theoretical considerations: this measurement defines a function g(x). 

[0100]The noise will thus be modelled by the composition of two operators: 

A transport operator, as in the first application example ; 

An amplitude modulation, that is the multiplication of the solution to the 
transport equation by function g(x). 

[01 01] If again the scheme given for the 1 st application example is used (and notably 
if the same hypotheses and . notations are used), the following information is 
introduced: 

a space Bj of noise-generating functions which will be the space of support 
functions j3j(t) in [tj min f tj max ] c [0,T] and that we extend by 0 in the whole interval [0J] ; 

noise modelling operator Tj: 

TjiPjMeBj^bjix.OeD 

solution to the transport equation 

Vb j {x,t)Jcj(x,t) = 0 in [ Xnia .x^ ]x [0,7] 
and meeting the initial conditions: 

f* y (*,i = 0) = 0 

[0102] In fact, the transport equation is solved by means of a numerical scheme. It is 
for example possible to use the conventional finite-difference centered scheme. 
Therefore discretization intervals Ax/ and Atj are selected (which are not necessarily 
equal to Ax and At, but which are assumed to be submultiples of these quantities, 
even if more general situations can be considered) and a grid is introduced whose 
nodes are the points of coordinates (x^ + i'AXj', n'Atj) with i' and n' e N. The 
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conventional finite-difference centered scheme explained hereunder allows, from the 
initial conditions, to calculate step by step the different values taken by function 
bj(x f t) at the various nodes of the grid (to simplify the notations, subscript j 
associated with the correlated noise considered has been omitted). 


1 \h n ' + l 4- h n ' + l - h n ' - f> n 'l 4- 


2A^ 

1 kn' + l , in' _ rn'+l _ .n'l _ f} 


2 Ax 1 r n+ 5 


[0103] In the above formula, quantity 04 "' represents the evaluation of quantity a (a is 
a generic term) at the point of coordinates (Xmm + i'AXj', n'Atj ). 

[0104] Here, discretization intervals Axj and Atj' have to be selected sufficiently small 
in relation to the wavelength because the numerical scheme is desired to give a 
precise approximation of the function solution to the transport equation (no 
dispersion is desired). 

[0105]The last stage of the noise modelling procedure selects the samples br n> that 
belong to the nodes common to the two grids and in multiplying them by g(iAx): the 
various components of the vector (of the data space) Tj(Pj) are obtained. 

[01 06] Specification of a norm or of a semi-norm in the data space 

[01 07] A first choice takes as the norm in the data space any discrete norm meant to 
be an approximation (via a quadrature formula) to the norm defined on a continuous 
basis by: 
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[0108]The semi-norm can also be used, which is a better choice as discussed 
below: 



where u is any vector of the data space (i and n are the indices representing the 
sample numbers respectively in space and in time). 

[01 09] Besides, other choices are also possible. 

[01 10] Specification, in each noise-generating functions space, of a norm for which 
each noise modelting operator establishes an isometric relation (or the 
approximation to an isometric relation) between the noise-generating functions 
space and the data space (provided with the semi-norm or defined) 

[01 11] Herein the procedure is specified that can be used to define the norm in the 
noise-generating functions space associated with each correlated noise. 

[01 12] A first choice provides each noise-generating functions space Bj with any 
discrete norm meant to be an approximation (via a quadrature formula) to the norm 
defined on a continuous basis by (here again, the subscript j has been omitted to 
simplify the notations): 
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T.max 


Pit) 


in which case, using the fact that bj(x,T)=0 for x e [XmirnXmaJ. linear operator Tj 
performs an isometry between Bj and D. However, since calculation of Tjpj is carried 
out numerically, the isometric character of operator Tj will not be perfect but only 
approximate, the approximation being all the better as discretization intervals Ax, At, 
AXj\ At/ are small. The fact that only approximately an isometry exists will lead to 
lower performances concerning the algorithmic implementation presented in 
paragraph 8. 

[01 13] A better choice provides space for each noise-generating function Bj with the 
semi-norm (here again, the subscript j has been omitted to simplify the notations): 


[01 14] In this case, operator Tj rigorously performs an isometry between Bj and D, at 
least if the solution to the discretization scheme of the transport equation is 
identically zero for n=N. 

[0115] Definition of the cost function 

[0116]The cost function is defined by the formula as follows: 


\\P\\b = + AarA* 



C(£m, 0 X9 fi 2 ,.../?, ) = (m 0 ) + F(5m) + fj.fi j - d\ 


Id 
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[01 17] It can be noted that, in the above function, F N L(m 0 ) is a constant vector. 

[01 18] See/ring the model and the noise-generating functions minimizing the cost 
function, this search being carried out by an algorithmic method taking advantage, 
explicitly or implicitly, of the isometry properties of the noise modelling operators 

[0119]The algorithmic method uses the specificity of Schur's complement, a 
specificity associated with the isometric character of operators Tj. This can be done 
for example by realizing that minimizing C(6m i p 1 , p 2 ,..., Pj) amounts to minimizing 

C(Sm) = CiSmJ^J.iSm^JjiSm)) 

where (p x (Sm) 9 p 2 {5m),...p J (6m)) minimizes C(8m,pi, p 2 ,..., pj) for given 5m. It can be 
noted that the Hessian associated with this quadratic form is Schur*s complement. 
C'(8m) can be minimized by means of any optimization method ; cost function C'(8m) 
being quadratic, the use of a conjugate gradient method is here particularly 
appropriate, 

[0120} Determination of (fi } (m) 9 fi 2 (m) i ...fi J (m) is very easy in view of the isometric 
character of operators T jf notably if choices (4) and (5) have been respectively 
selected to define || \\ D and || (J . If there is only one type of noise (J=1 ), which is the 

case for our illustration, and if the isometry is perfect (which is also the case thanks 
to the choices made), determination of P x (m) simply requires evaluation of the 
gradient of the quadratic form. If there is a superposition of several noises (J>1), 
various algorithms can be considered to take advantage of the isometric character of 
operators Tj as explained for the first application example. 
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[01 21] Figure 14B shows the results obtained (seismic response of the solution 
model) by implementing the method when function g(x) is selected in a very simple 
form (affine function). The improvement in relation to the conventional inversion 
result can be observed (Fig.14A). 

[01 22] Implementation examples where the physical parameter whose subsurface 
distribution is modelled is the acoustic impedance have been described. It is clear 
that the method, in its most general definition, can be applied for seeking the 
distribution of physical quantities affected by correlated noises in any heterogeneous 
medium. 
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